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Abstract 

Background: The interest in cell membrane has grown drastically for their important 
role as controllers of biological functions in health and illness. In fact most important 
physiological processes are intimately related to the transport ability of the membrane, 
such as cell adhesion, cell signaling and immune defense. Furthermore, ion migration is 
connected with life-threatening pathologies such as metastases and atherosclerosis. 
Consequently, a large amount of research is consecrated to this topic. To better 
understand cell membranes, more accurate models of ionic flux are required and also 
their computational simulations. 

Results: This paper is presenting the numerical simulation of a more general system 
modelling ion migration through biological membranes. The model includes both the 
effects of biochemical reaction between ions and fixed charges. The model is a 
nonlinear coupled system. In the first we describe the mathematical model. To realize 
the numerical simulation of our model, we proceed by a finite element discretisation 
and then by choosing an appropriate resolution algorithm to the nonlinearities. 

Conclusions: We give numerical simulations obtained for different popular models of 
enzymatic reaction which were compared to those obtained in literature on systems of 
ordinary differential equations. The results obtained show a complete agreement 
between the two modellings. Furthermore, various numerical experiments are 
presented to confirm the accuracy, efficiency and stability of the proposed method. In 
particular, we show that the scheme is unconditionally stable and second-order 
accurate in space. 

Keywords: Reaction-diffusion system, Electromigration, Nonlinear coupled system, 
Finite element method, Nernst-Planck equations, Numerical analysis, Enzyme kinetics, 
Substrate suicide, Cooperative phenomena, Computational simulation 



Background 

Cell membrane is the biological membrane separating the intracellular environment from 
the extracellular one. The cell membrane surrounds all cells and it's selectively permeable, 
permitting the free passage of some substances and restricting the passage of others, thus 
controlling the flux of substances in and out of the cell. All diseases are problems of reg- 
ulating the passage of materials at the level of the cell. Consequently, to understand the 
cause of a disease, we need to understand the alterations that take place at the cellular 
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level. Thanks to mathematical modeling, cell biology phenomena may be expressed by 
ordinary differential equations or systems of partial differential equations. 

An important class of models in cell biology, is the class modeling ion transport through 
biological membranes. This transport phenomena occurs in most living cells and some 
biochemical processes. The first models in the literature included one ordinary differ- 
ential equation for each ion concentration. All this models were based on the implicit 
assumptions that chemical concentrations are uniform in space. This assumption is rea- 
sonable when the region of space where the reaction occurs is confined and very small. 
Also this models assumed that the electric field is constant inside the membrane. 

Despite that the constant electrical field assumption has the advantage of leading to 
a simple mathematical analysis, all cells maintain a difference in electrical charge across 
their membrane. This difference in charge give rise to a voltage difference, or electrical 
potential. Furthermore, there are numerous situations in which chemical concentrations 
are nonuniform in space. In this sense, we need to establish and compute more accurate 
mathematical models of ions electromigration through biological membranes. 

In this paper, we present the numerical simulation of a more realistic model of ions 
electromigration through biological membranes. This model is more general than those 
in literature of membrane transport as it extends them in four topics: 1) it's a multidimen- 
sional model, 2) it doesn't rely on the constant electrical field assumption, 3) it consider 
both the temporal and spatial dependence, 4) it includes different reaction kinetics terms. 

Introduction 

In this paper we consider a class of models of ions migration through biological mem- 
branes. Such migrations exist for most living cells and some biochemical processes. The 
motion of ions is supposed due to diffusion and to the effect of the electrical field. 
Furthermore ions can undergo reactions. So the ions concentrations satisfy the Nernst- 
Planck equations, including a kinetic reaction terms and the potential is given by Poisson 
equation. The model is given by: 



3Q 
dt 



- diACi - midiv(CiV4>) = Fi(C\, C Ns ) in Q T , for i=l,...,Ns 



- eA0 = y^ZjCj -f 

, 3Q 30 
di— +rntCi— = 0 
dv dv 

<p(t,x) = 0 
Q(0,x) = C ii0 (x) 

0(0,*) = 00 (X) 



in Qt 

in for i = 1, . . . , Ns 

on £2 
on Q 



(1) 



where Qt = ]0, T[ x £2, J] r = ]0, T[ x T > 0; £2 is an open regular set of K 2 which 
represent the biological cell and dQ represent the cell membrane. 

For each i, Q is the concentration of the i species which has diffusion coefficient di, 
mobility m; and valency z;. 0 is the electrical potential,/ is the fixed charges concentration 
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and the F; are reaction terms. We suppose that F/ depends continuously on the Cj 's and 
4>, and that / is a bounded function. We suppose that di is a positive constant for each i. 



Further information about the modelling of the problem and its mathematical analysis, 
can be found in [1]. 

In the biophysical literature, the early works on these models were interested in the 
stationary case of passive migration (i.e., without reaction); see [2,3]. In all these works 
two popular simplifications were considered, namely the Goldman hypothesis where 
the electrical field is supposed to be constant inside the membrane and the electroneu- 
tral hypothesis where the neutrality at each point of the membrane is assumed; see 
for example [3]. Mcgillivray [4] recognized that these models are the limit of the full 
equations when the ratio *J~e = j of the Debye length to the membrane thickness goes to, 
respectively, infinity or zero. Usually enzymes are held to biological membranes and ions 
undergo biochemical reactions when crossing the membrane. Valleton [5] did a general 
biophysical study of coupling of electromigration diffusion with biochemical reactions. 

In this paper we present a numerical simulation of such systems, for a large class of 
reaction kinetics, including the usual biochemical kinetics as the Michaelis-Menton one 
(a mathematical analysis of the one dimensional and stationary case was done by [6] 
then we did the mathematical analysis of the multidimensional unsteady case [1]). This 
article is organized in the following way. The next section is devoted to finite element 
discretisation of the mathematical model. Then, we present applications, results and 
numerical experiments showing the accuracy, efficiency and stability of the proposed 
method. Finally, conclusions are drawn in the last section. 

Finite element discretisation 

In order to show the numerical formulation of the problem, let V = L 2 (0, F;// 1 (S2)) be 
the space of approximate solutions and W = H (Q) be the space of tests functions. Let 
W h be a finite element space of Lagrange PI included in W and V h = L 2 (0, T; W h ) be 
the finite dimensional subspace of V. Now we introduce the function Zp, = (Zi : h)i<i<N s 
defined by 



The Faedo-Galerkin formulation for the problem is given by, finding Zy, € V/, for i = 
1,...,N S and <f>h € Vf, such that (/>/,= 0 in 3Q : 

• for every Wh € Wh a.e. t e] 0, T[ and for 1 < i < Ns, 



for all i = 1, . . . ,Ns Q,o e L 2 (Q.) and satisfy Q,o > 0 



(2) 




Moreover, we consider 




ftfn <H,h z i,h w h + f Q dtq ith VZ ith Vw h = f a Fj(Z h )w h 
Z i>h (0) = ?4*(0)Cjb, 



• for all Vf, € Wf, such that V/, = 0 in 3£2 and a.e. t e] 0, T[ 

s fn V(/>* vv fc = f a (E ZjqjM z j,h -f)n 



(3) 



/•=1 



4> h (0) = $ onQ 
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where Fi(X\, . . . ,X Ns ) = Fi(q 1>n X lt . . .,qN s ,h x N s ) and C* 0 , $ are the projections of 
Q,o, <j>o on W n . 

Let (y/)i<j< m the mesh nodes, (®j)i<j< m the canonical basis of W n , we consider the 
following two sets of index 

Yo = {/e [l:m], yj e } 
Y=[l:m]-Y 0 

m 

We represent the solutions as Z is! fo(t,x) = J2 a Z isi j(t)®j(x) for i s = 1,...,N S and 

/'=! 

m 

0 fc (t,«) = E«y(*)*y(«) with a/ (f) = 0 for; € Y 0 . Then we set fe. = fejt) = 
/=i 

(«z is ,/(0)i</<m for z's = 1, . . . ,A4 and ^ = ^ (f) = (<*/(*)) ls/;£> „. 

Now let's consider an uniform subdivision of [0, T], we define a time step dt = jj, for 
N > 1. We pose then: 

= «<sft, 0 < n < N 

Let's note fz« the approximation of fz is (£„) and ^» the approximation of ^(t„), then we 
used an implicit scheme for the discretization of the time derivative. By the method of 
finite elements see [7], we arrive at the following formulation of the problem: Find the 

vectors nodal concentrations £ 7 «+i = (a 7 n+i .) for every 1 < i s < N s and nodal 

V 's •'/l<J<m 

potential = (a" +1 ) with a" +l = 0 for / e Y 0 such that 

v \ I ) \<j<m 1 



<) 

with af +1 = 0 foi 

•j<m 

_ 1 A . t . A . . t 

dt 



£ z o =fe is (0), f^ = f#(0) 

Let T n the mesh generation of Q. containing nel finite elements. For a finite element 
€ T n , let be T e * = fe} where k\, ki, k% are the numbers of degrees of free- 

dom of and N%j are interpolation functions. We have <&i/e k = Nj£, ^>j/e k = where 
s, / € {1, 2, 3} if and only if (i, /) e T e * x 7 e *. 

Where = (4")^^ = (ff)^ ^ W ( V 1 ' ' " ' ' %0 = 

(5; s '" +1 ) for 4 = 1, . . . ,A/ S , Q = (qy)i<y< m , M isi( = (n/f) ,fi = 

\ / l<i<m *^ \ ' / l<i,j<m 

(Pi)i<i<m (XyxY means the extracted matrix from X by keeping lines and columns with 

numbers belonging to 30 



nel „ / m \ we/ 

r = E/ «p ^E«f% = E4' H " 



where 



h,n,e k 
a ij 



r h,n 
i) 



f exp (^E a . 

J <>k \ U 's p= l 



IK if (i /) e ^ x r* 



0 ifior/gr* 

tiel yyi tiel 

2> / ^<^E«r*^%™<* = EC 
*=i ^ !s /=i *=i 
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where 



h,n,e k 
i) 



dis L exp =2* E « VATf* V< if (i, ;) € ^ x r* 



E / ^ 



if / or / ^ T ek 

n 

ex p ( x^ a " +1$iie * ) ■ ■ ■ • 



we/ 

— W s, " +1,e * 



/=1 



*=1 



where 



0 



where 



Me/ . ne/ 

E/ v % v *«k* = E«# 



0 if z or / £ T e * 



r=E/ 



-p ( ^E-r*'* ) = E<" 



<4 



i,e k 



where 



m 



k,n,e k 



4 exp ^ £ «« A* A#V* if (i, j) e T°« x r* 



if ior; £ r e * 



where 



f /a£* ifze T e * 
J e<r 



o if ;g r e * 

Finally, finding £ 7 „+i = \a 7 n+i •) for every 1 < z' s < A/ s and f^+i 

is V k '' / Ki<m 



(a" +1 )i</< m with = 0 for ;' € Y 0 such that 



H z o = fej,(0), ? 0 o = ^(O) 



We have a nonlinear term due to Sj s ,^„ +1 (£ z k+i, • • • > §z«+0 > we nave dealt with according 
to the model and thus to the expression of +1 . 
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Results and discussion 

In this section we present three numerical applications of ions electro-migration through 
biological membranes. The models of the basic enzyme reaction, the suicide substrate 
reaction and the cooperative reaction, are numerically simulated. 

Result 1 : Enzyme kinetics (basic enzyme reaction) 

To understand where some of the more complicated reaction schemes come from, we 
consider a reaction that is catalyzed by an enzyme. Enzymes act as remarkably efficient 
catalysts (generally proteins), by accelerating the conversion of some other molecules 
called substrates into products via lowering the free energy of activation of the reaction, 
but they themselves remain unchanged by the reaction. Thus, they are important in the 
regulation of biological processes, for example as activators or inhibitors. One of the most 
basic enzymatic reactions, first suggested by Michaelis and Menten [8], implies a sub- 
strate S reacting with an enzyme E to form a complex SE which in turn is converted into 
a product P and the enzyme. This is represented schematically by 

S + E^SE^P + E (4) 
k-i 

here k\, k-\ and k 2 are constant parameters associated with the rates of reaction. We 
denote the concentrations of the reactants by 

Ci = [E] , C 2 = [S] , C 3 = [P] , C 4 = [ES] . 

Then the law of Mass Action applied to (4) leads to one equation for each reactant 
and hence a system of nonlinear equations. The usual approach to these equations is to 
assume that the initial stage of the complex C4, formation is very fast after which it is 
essentially at equilibrium, then we get C4 in terms of C 2 , 

„ Ci,oC 2 , k-\ + k 2 

C4 = - - - , k M = : 

C 2 + K M k\ 

The basic enzyme reaction model becomes then 
9Ci 



d\ ACi — m\div(C\Sl '<j>) = 0 in Qj 

— 1 - d 2 AC 2 - m 2 div(C 2 V<l ) ) = ~cI+Km in 

^ - d 3 AC 3 - m 3 div(C 3 V<l>) = in Q T 

3 

-eA0 = J>;Q -/ in Q T 

dCi '~ 30 ^ 

di- YmiCi— =0 in 2_,t» for i = 1,2,3. 

dv 3v 

(j)(t,x) = 0 in J^t 

Ci(0,x) = Ci,o, C 2 (0,x) = C 2 fi, C 3 (0,x) = 0 on Q. 

(j)(0,x) = <po(x) on £2 

(5) 

Algorithm of resolution 

Before stating the resolution algorithm, we introduce the function Zf, = (Zi,h) 1<i<3 
defined by 



Zi,h = C uh exp (^f(<l>h)j for i = 1, 2, 3 
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Moreover, we consider 



p uh = exp \"^<t>h) and^y, = — 
V di J Pi, h 



We used the following algorithm to calculate </>/, and Zy, then we calculate Cy by using 
the reverse relation: 

(—mi \ 

• Initialize: for i = 1, 2, 3, 

Z lh = Ci i0 (0,x)p ith (0,x), 

0 (—mi 
1i,h = ex P ( —r-4>h(0,x) 



d. 



• Loop over n 
At step n : 



Calculate </>f +1 solution of 



• Calculate ^J 1 , ^J 1 by: 



• Calculate Z"^ 1 solution of 



r a n+l 7 n+1 - a" 7" r 

/ — x : ~w h + d x \ (flX VZJ Vwj = 0 

Jo. at Jq ' 



• Calculate Z^ 1 solution of: 



z „+u+i_ z „+u 



Loop over k untill 
where eps is the stopping criterion. 



Li (fi) 

n+l 7 n+l,k+l „n 7 n „ „ z._ _ „B+1 7 n+l,k+l 

n+l 7 8+U ^ 



• Calculate solution of 

<l3,h Z 3,h -l3,h Z 3,h , j f „+i„^„+i„ f k 2 C h0q 2 ,h Z 2,h 



f %,h Z 3,h %3,h 3,h , . /" »+l n7 »+l n /" 



Numerical results 

Here we present changes in substrate, product and enzyme concentrations. The cell is 
represented by an ellipse with semi-major axis a=2 and semi-minor axis b=l. The dif- 
fusion coefficient of the ions are d\ = 10 -3 m 2 .s -1 , d 2 = 2.10 _3 m 2 .s _1 and d-$ = 
5.10~ 3 w 2 .s _1 . The constants of reaction are ku = 9.10 -5 M, k 2 = 1,4.10 4 s -1 
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and k\ = k cat = j^- = 1, 55.10 8 M~ l .s~ 1 . The charge number of the ions are z\ = 1, Z% = 
0 and Z3 = 1. The electric charge density is/ = 0.1C. The initial conditions are Ci,o = 
lfiM, C2,o = 800 fiM, Czfl = 0 and 4>o = —80 mV; the stopping criterion is eps = 10 -4 . 
The time step of the simulation is dt = 10 _3 s. 

Shown in Figure 1 is the spatial distribution of substrate concentration through the cell 
at both initial and final time (t=0 and T=200 ms). 

Figure 2 presents the spatial distribution of the product concentration through the cell 
at both initial and final time (t=0 and T=200 ms). 

In Figure 3, one sees the time evolution of the substrate and the product concentrations 
at the center of the cell. It can be seen that the substrate decrease curve is the mirror 
image of the product appearance curve. By observing the early times, it's obvious that 
the substrate loss and product appearance change speedily with time but as time goes on 
these rates diminish, to reach zero when all the substrate has been converted to product 
by the enzyme. 

The Figure 4 shows, as predicted, the enzyme concentration remaining constant over 
time. 

Two simplifications of these equations have been quite popular in the literature while 
computing membrane reactions, firstly the Goldman hypothesis where the electrical field 
is supposed to be constant inside the membrane and secondly considering a system of 
ordinary equations depending only on time and not the space. The added value of this 
work is not considering all of those simplifications which leads to a more realistic model 
and more accurate numerical results. Moreover, the results obtained are in agreement 
with the experimental results found in the literature [9]. 

Result 2 : Suicide substrate kinetics 

An enzyme system of major experimental concern; see [10,11], is the mechanism-based 
inhibitor, or suicide substrate system, represented by Walsh et al. [12], 

(6) 

Y k -\E t 

where E, S and P stand for enzyme, substrate, and product, respectively; X and Y, enzyme- 
substrate intermediates; £,-, inactivated enzyme; and the k 's are positive rate constants. 

In this system, Y has a choice of one of two pathways, namely, to E + P with rate 1<3 
or to Ei with rate fcj.. The ratio of these rates, /C3//C4, is called the partition ratio and is 





t=0 seconds — '«™ t=0.2 seconds MHWll 

Figure 1 The spatial distribution of the substrate concentration through the cell at initial time t=0 
and final time T=0.2 s. 
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t=0 seconds 




t=0.2 seconds 



Figure 2 The spatial distribution of the product concentration through the cell at initial time t=0 and 
final time T=0.2 s. 



denoted by r. Each of these pathways are supposed to be irreversible over the timescale 
of the reaction see [13]. S is known as a suicide substrate because it binds to the active 
site of an enzyme — like a substrate — but the enzyme converts it into an inhibitor which 
irreversibly inactivates the enzyme. Thereby, the enzyme 'commits suicide'. In this way, 
a suicide substrate can specifically target an enzyme for inactivation. Furthermore, sui- 
cide substrates are particularly useful in drug administration, as they are not noxious in 
their common form and only the designated enzyme can convert them to their inhibitor 
form. For example, suicide substrates have been subject of investigation for use in the 
treatment of depression (monoamine oxidase inhibitors, Seiler at al. [10]), epilepsy (brain 
GABA transaminase inhibitors, Walsh [11]), and some tumors (ornithine decarboxylase 
inhibitors, Seiler et al. [10]). Suicide substrate kinetics have been studied by Waley [13] 
and by Tatsunami et al. [14], who had interest in the factor which determined whether 
the substrate was exhausted before all the enzyme was inactivated. Waley proposed it was 
rfi, where fi is the ratio of the initial concentration of enzyme to that of substrate, namely, 
en/so- Tatsunami et al., on the other hand, found the determining factor to be (1 + r)(i. 
When (1 + r)fi > 1 the substrate is exhausted, while for (1 + r)/i < 1, all the enzyme is 
inactivated. When (1 + r)ji = 1, both occur. The interest is when eo/so is not small, which 




"substrate.dat" 
"product.dat" 



0 20 40 60 80 100 120 140 160 180 200 

Time (ms) 

Figure 3 Evolution in time of substrate and product concentrations at the center of the cell. 
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80 100 120 
Time (ms) 



1 80 200 



Figure 4 Temporal evolution of the enzyme concentration at the center of the cell. 



was in effect assumed since both Waley [13] and Tatsunami et al. [14] used a quasi-steady 
state approximation. The validity decreases for increasing values of eo/so- We denote the 
concentrations of the reactants by 

Ci = [E] , C 2 = [S] , C 3 = [X] , C 4 = [Y] , C 5 =[Ei] , C 6 =[P] 

The law of mass action applied to (6) leads to one equation for each reactant and hence a 
system of nonlinear equations. We obtain the following system 

3Ci 
9C?2 _ 

i; 

3C?4 

it 

dt 

-sA<p 
dCi 

di h miCi — =0 in V r , for i = 1,2,3, ... ,6. 

(p(t,x) = 0 in J^t 

Ci(0,x) = eo, C2(0, x) = so, C,(0,*) = 0 on Q,, for i = 3, . . . , 6. 

<p (0, x) = <po (x) on £2 

(7) 

Algorithm of resolution 

Before stating the resolution algorithm, we introduce the function Z/, = (Z,' j / ! ) 1<i<6 
defined by 

ffli 

Z it h = C uh exp(— (<£/,)) for 2 = 1, .... 6 



rflAQ — midiv{C\V<p) 


= -hc x c 2 + 


/c_iC 3 + /c 3 C 4 


in Qj- 


d 2 ACi — rri2div{C2^ 4>) 


= -kidC 2 + 


/c-iC 3 


in Qr 


<f 3 AC 3 — msdiviCsV (p) 


= k x C x C 2 - (k 


-1 + /c 2 ) C 3 


in Qj™ 


d^ACn — m^diviCiVc/)) 


= k 2 C 3 - (/c 3 4 


- /C4) C4 


in Qt 


d5tS.Cs — msdiv(Cs\ '(p) 


= /C4C4 




in Qj™ 


deAC6 — m(ydiv{C(^ <p) 


= /C3C4 




in Qj- 


= E«iQ-/ 
i=i 






in Qt 


+ m,<~/ =0 
ov 




in for i 


= 1,2,3, 
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Moreover, we consider 

Pi,h = exp(^0 A ) andg-y, = — 
«i Pi,h 

We used the following algorithm to calculate (/>/, and Zy, then we calculate Cy, by using 
the reverse relation: 

Q,h = ex P(—} —4>h)Z i? h 
di 

• Initialize for i= 1, ... ,6 



z lh = Ci,o(0,x)p i>h (0,x), 

n / —mi \ 

%h = ex P \—jr<l>h(0-x) J 

• Loop over n 
At step n : 

• Calculate (/>^ +1 solution of 

• Calculate a n+l a n+l a n+l a n+l a n+1 a n+1 hv 
^aicuiaie q yh , q 2 h , q 3h , q 4h , q 5 h , q 6 h oy. 

• Calculate Z"+ h \ Z*+\ Z"+ h \ Z"+ h \ Z*+\ Z^ 1 solutions of: 

- initialize Z"^ 1 ' 0 = Z" iM for i = 1, .... 6, 

i.i-l 



Loop over Jt untill £ \^t* +X ~ Zf?* 



L 2 (Q) 



r a n+l 7 n+1 - k+1 - a n 7 n r 

L it wh + di L vz ^ h VWh 

r (7 «+ 1 7"+ 1 ^+ 1 - a n Z n r 

A * w * + d2 L *** vz ^ VWh 



k+1 7 n+l,k+l n -yn 



f ? 3,/i Z 3,ft ' g3,A Z 3,fe „, , /" «+l „„+U+l 

y Q — ^ — w * + * L % » vz ^ Vwh 

y Q * «* + * y s ^ vz 4> , 

= jf (AmJJ 1 ^ 1 - (As, + fc) <f Z^* +1 ) Wh 
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r a n+Y Z n ^' k ^ - a" Z n r r 

J Q J t Wh+ds J^q Sih VZ 5h Vw h = J Ua^ h Z %h w h 

r a n + l z n+hk+l - a" Z n r r 

J J t W h + d 6 J q VZ 6* V w h = J q *3«S ^ W ft 

Numerical results 

Here we present changes in substrate, product, enzyme, inactivated enzyme and 
the intermediate concentrations (Xandi^). The cell is represented by an ellipse 
with semi-major axis a=2 and semi-minor axis b=l. The diffusion coefficients 
of the ions are d\ = 10" 3 m 2 .s"\ d 2 = 2.10" 3 m 2 .s" 1 , d 3 = 5.10" 3 w 2 .5 _1 , J = 
10 _3 w 2 .5 _1 , ds = 2.10 _3 w 2 .s _1 , de = 4.10 _6 w 2 .s _1 , the reaction parameters are ki = 
2 s _1 , /c_i = 4 s _1 , ^2 = 12 /C3 = 10 s _1 and /C4 = 2 s _1 . The charge number of the 
ions are Z\ =\, Z% = 0, Z3 = 1, Z4 = 1, Z5 = 1 andz6 = 0. The electric charge density is 
/ = 0.1 C. The initial concentrations are en = 0.5/j.M andsn = 0.5 fiM; and <po = — 80 mV. 
The time step of the simulation is dt = 10 _2 s. The data employed for the reaction 
parameters and initial concentrations were taken from Burke et al. [15]. 

Figure 5 plots the changes in the concentration distribution of substrate from initial 
time to final time. One can see that in final time the substrate was totally exhausted. This 
complete consumption of the substrate is in agreement with the prediction of Tatsunami 
et al. [14] as (1 + r)/x = 6 > 1. 

In Figure 6, one sees the evolution of substrate concentration over time at the center of 
the specimen. 

Figure 7 shows the evolution of inactivated enzyme concentration over time at the 
center of the specimen. 

Figure 8 shows the numerical solutions for intermediate concentrations of X and Y. 

In Figure 9, is represented the graphic of the evolution of substrate and product at the 
center of the cell, comparing that result with Figure 3 (Michaelis and Menten model), here 
the two plots are asymmetric which is logical as we know that an amount of the substrate 




Figure 5 Substrate concentration from initial time to final time. 
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Figure 6 Evolution of substrate concentration overtime at the center of the cell. 







instead of being converted to product, is forming the inactivated enzyme (inactivating the 
enzyme). 

To illustrate, Figure 10 shows the decrease in the enzyme concentration unlike the 
Michaelis Menten model (Figure 4); however, as the intermediate enzymes X and Y, van- 
ish in few milliseconds, we see the loss in enzyme compensated by the production of the 
inactivated enzyme: the enzyme commits suicide. 

To highlight the accuracy of these results, we compared them first with the numer- 
ical solutions and the approximate asymptotic solutions obtained both by Burke et al. 
[15], they considered a system of ordinary differential equations depending on time as 
they neglected the spatial aspect of the biochemical reaction, and they supposed that 




Time (1/100) 

Figure 7 Evolution of inactivated enzyme concentration over time at the center of the cell. 
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Figure 8 Evolution of the intermediate concentrations X and V over time at the center of the cell. 



the electrical potential inside the membrane remains constant. For the numerical solu- 
tions Burke et al. [15] solved the system numerically, but for the approximate asymptotic 
solutions they non-dimensionalise the same system, and used asymptotic methods and a 
method detailed in Kevorkian and Cole [16]. Finally we compared our results with pre- 
vious approximate methods of Tatsunami et al. [14] and Waley [17] which were based 
on a pseudo-steady state hypothesis. This comparison shows that the results described 
here are valid numerical solutions for the kinetics of suicide substrate system. The 
solution for the substrate and inactivated enzyme are more accurate than those of pre- 
vious approximations [14,17], especially in small time, which is by definition ignored 
by any pseudo-steady-state approximate method. Furthermore, the method presented 




Time (1/100) 

Figure 9 Evolution of substrate and product concentrations over time. 
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Figure 1 0 Evolution of enzyme and inactivated enzyme over time at the center of the cell. 



here is specially useful in estimating the intermediate (X and Y) concentrations besides 
incorporating the spatial and the electro-migration aspects of the phenomena. 

Result 3 : Cooperative phenomena 

An enzyme reaction is said to be cooperative if a single enzyme molecule, after binding a 
substrate molecule at one site can then bind another substrate molecule at another site. 
Such phenomena are quite common in living organisms. Another interesting cooperative 
reaction is when an enzyme with several binding sites is such that the binding of one 
substrate molecule at one site can affect the activity of binding other substrate molecules 
at another site. This indirect interaction between distinct and specific binding sites is 
called allostery, and an enzyme displaying it, an allosteric enzyme. When a substrate that 
binds at one site increases the binding activity at another site then the substrate is called 
an activator, otherwise (if it decreases the activity) it's called an inhibitor. 

As an example of cooperative phenomenon we consider the case when an enzyme has 
two binding sites. A model for this consists of an enzyme molecule E which binds a sub- 
strate molecule S to form a single bound substrate-enzyme complex X. This complex X 
not only breaks down to form a product P and the enzyme E again, it can also combine 
with another substrate molecule to form a dual bound substrate-enzyme complex Y. This 
Y complex breaks down to form the product P and the single bound complex X. The 
reaction mechanism is represented schematically by 



s- 


f-£ 


!± 


X 




S -\ 


-X 


*~3 


Y 





(8) 

X + P 



Here k's are rate constants. We denote the concentrations of the reactants by 



d = [E] , C 2 = [S] , C 3 = [X] , C 4 = [Y] , C 5 = [P] . 
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Then the law of Mass Action applied to (8) leads to one equation for each reactant and 
hence the system of nonlinear reaction. We have then 

3 C\ 

— d\ ACi - midiv(CiV<p) = -k\C 2 C\ + (k-i + k 2 ) C 3 inQj- 

i£ 2 



- d 2 AC 2 - m 2 div(C 2 V <j>) = -k\C 2 C\ + (k-\ - k 3 C 2 ) C 3 + £_ 3 C 4 in Q T 

i6 3 



d 3 AC 3 - m 3 div( C 3 V0) = hC 2 Ci - (k-i + k 2 + k 3 C 2 ) C 3 + (k- 3 + £ 4 )C 4 in Q T 
— d 4 AC 4 — m^diviCn^i '(p) = k 3 C 2 C 3 — (k- 3 + k 4 )C 4 in Qj 



9(^4 

— d s AC 5 - m 5 div(CsV(p) = k 2 C 3 + kiC 4 in Q T 

at 

5 

-sA<p = £ ZjQ -/ in Q T 

dQ d<j> _ „ 

flj h WiQ — = 0 in 2^t, for j = 1, . . . , 5 

4>(t,x) = 0 in J] r 

Ci(0,^) = eo, C2(0,x) = so, Q(0,x) = 0 on f2, for j = 3, 4, 5. 

0(0, x) = 4>q(x) on £2 

(9) 



Algorithm of resolution: 

Before stating the resolution algorithm, we introduce the function Zj, = {Zi,h) 1<i<5 
defined by 



Zi,h = C uh exp (jj;(<i>h)j for / = 1, . . . , 5 



Moreover, we consider 

Pi,h = exp ("j-<Ph) and q uh = — 

We used the following algorithm to calculate 0/, and Z ;> /, then we calculate Q,/, by using 
the reverse relation: 

(-»!; \ 
4>h I z, )A 

• Initialize for z = 1, . . . , 5 

z y< = Ci,o(0,x)pij,(0,x), 

ilh = ex P y-^'Phio.x) 

• Loop over n 
At step n : 

• Calculate 4>^ +1 solution of 

5 



. Calculate^', by: 
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• Calculate Z£+\ Z£+\ Z n + h l , Z n + h l , Z n + h l solutions of: 
- initialize Z*+ 1,0 = Z n ih for i = 1, .... 5, 



- Loop over k untill J2 

i=l 



Z' 



n+l,k+ 1 



£2 (£2) 



«+l 7 n+l,k+\ 



a" 7 n r 

./£2 ' 



Jn 



y Q — ^ — w " + j2 y Q *** vz ^ vw " 



h+1 7 n+l,k+l 



Jn 

/ 



a" 7" 
^.h 3,h 



dt 



w h +d 3 



f < 

Jn 



v2 jH-U+l Vwfc 



.n+1 7 n+l,k+l 



= J a (^?2,r^ u+l <r^ u+1 -(*-i+*2+* >c 1 



/• a" +1 Z" +1 ' k+1 - a" Z" r 

y Q — ^ — w " + j4 y Q *** vz ^ vw " 

= / q (* S ^ 1 ^ 1 ^ 1 ^ 1 -(*-3+*^ 1 ^ 1 ) 



/ 

Jn 

/ 



h+1 7 n+l,k+l 
5,h 



a" • 7' 



n n 7 n 
^5,h 5,h 



Wh + ds 



Jn 



dt ~" • " J J Q ^ h 5 * 



„ +lvz „+ u+ i v ^ 



Numerical results 

Here, we present the evolution of product and enzyme concentrations over time in both 
positive and negative cooperativity. We used the same constant parameters as the previ- 
ous example for the diffusion coefficients of the species and the same initial conditions. 
The cell is represented by an ellipse with semi-major axis a=2 and semi-minor axis b=l. 
To ensure cooperativity, the positive rate constants are chosen by the following reason- 
ing: Suppose that the binding of the first substrate molecule is slow, but that with one 
site bound, binding of the second is fast (this is large cooperativity). This can be mod- 
eled by letting £3 -> 00 and k\ -> 0 while keeping k\k 3 constant, in which case K2 -*■ 0 
and Ki 00 while K2K1 is constant ( K\ and/<2 were introduced as they appear in the 
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Figure 1 1 Evolution of the product concentration in both positive and negative cooperativity. 



expression of velocity reaction, this is discussed in details in the book by Keener and 
Sneyd [18]), where 



Ki 
K 2 



k-i + k 2 

h 

/c 4 + k-3 



An enzyme can also exhibit negative cooperativity, in which the binding of the first 
substrate molecule decreases the rate of subsequent binding. This can be modeled by 
decreasing £3. We used for positive cooperativity K\ = 1000, K2 = 0.001 and K\ = 0.5, 
K2 = 100 for negative cooperativity (this values were taken from [18]). 
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Time (1/100) 



Figure 1 2 Evolution of the enzyme concentration in both positive and negative cooperativity. 
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Table 1 Convergence results for the basic enzyme reaction 

Mesh size h\= 0.3 fe 2 = 0.1 h$= 0.05 

Error 268.10~ 7 72.10~ 7 45. 10~ 7 



Table 2 Convergence results for the suicide substrate reaction 

Mesh size h\= 0.3 fe 2 = 0.1 hi= 0.05 

Error 93.1 0" 3 22.1 0" 3 56.1 0" 5 



Table 3 Convergence results for the cooperative reaction 



Mesh size 


/«!= 0.3 


h 2 = 0.1 


h 3 = 0.05 


Error 


49.10" 4 


28.1 0" 4 


93.1 0" 5 



In Figure 11, one sees that in positive cooperative reaction, the product concentration 
is characterized by an "S-shaped" sigmoidal curve, which is different from other enzyme 
reaction that exhibits a curve that tends to be hyperbolic. This results from cooperative 
effects; in which the enzyme can bind more than one substrate molecule, but the binding 
of one substrate molecule affects the binding of subsequent one. 

In Figure 12, we plot the evolution of the enzyme concentration for both extreme 
positive cooperativity and negative cooperativity. 

The convergence test 

The rate of convergence of the scheme is difficult to prove analytically. However, numer- 
ical experimentation suggests that the scheme is second-order accurate in space. A 
quantitative estimate of the convergence error was obtained by performing a number 
of simulations for the same initial condition on a set of increasingly finer space meshes 
and time steps. The initial conditions are constants. Let the mesh generation of Q, 
and h(T^) = max{diam{ei ( )\ei <: e we take h = 0.3, h = 0.1 and h = 0.05. For 
each mesh we integrate to time T with dt = Note that as we refine the space 
step we also refine the time step. The error of the numerical solution was defined as 

E (h) = dt x E £ max ||z^ +1 - Z*f | 

;=1«=0 k " llZ/(Q) 

The convergence of the basic enzyme reaction 

In Table 1 is presented the error of convergence for different mesh sizes in the case of 
basic enzyme reaction. 




(3) Time step dt= 0.00025 ( b) Time step dt=0.0005 (c) Time step dt=0.001 

Figure 1 3 Snapshots of the product concentration in the basic reaction at T=0.5 with three different 
time steps, shown in subfigures (a), (b) and (c). The time steps are given below each subfigure. 
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(a)Timestepdt=0.0025 (b) Time step dt=0.005 (c) Time step dt= 0.01 

Figure 14 Snapshots of the product concentration in the suicide substrate reaction at T=4 with three 
different time steps, shown in subfigures (a), (b) and (c). The time steps are given below each subfigure. 



The convergence of the suicide substrate reaction 

In Table 2 is presented the error of convergence for different mesh sizes in the case of 
suicide substrate reaction. 

The convergence of the cooperative reaction 

In Table 3 is presented the error of convergence for different mesh sizes in the case of 
cooperative reaction. 

Stability and accuracy tests 

Now, let us give some information about the numerical stability of our algorithms. We 
perform a numerical experiment with different time step dt, y and These results sug- 
gest that the scheme is indeed unconditionally stable as the solutions are quasi the same 
for different time steps. To illustrate, we chose to represent the product concentration. 

Stability of the basic enzyme reaction 

In Figure 13, We display snapshots of the product concentration at time T=0.5 with three 
different time steps dt=0.00025, dt=0.0005 and dt=0.001. We can see that the results are 
quasi the same at the final time T. 

Stability of the suicide substrate reaction 

In Figure 14, we display snapshots of the product concentration at time T=4 with three 
different time steps dt=0.0025, dt=0.005 and dt=0.01. We can see that the results are quasi 
the same at the final time T. 

Stability of the cooperative reaction 

Figure 15 shows snapshots of the product concentration at time T=4 with three different 
time steps dt=0.0025, dt=0.005 and dt=0.01. We can see that the results are quasi the same 
at the final time T. 




(a) Time step dt=0.0025 (b) Time stepdt=0.005 (c) Time step dt=0.01 

Figure 15 Snapshots of the product concentration in the cooperative reaction at T=4 with three 
different time steps, shown in subfigures (a), (b) and (c). The time steps are given below each subfigure. 
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Conclusion 

In this paper, a new model simulating ions electro-migration through biological mem- 
branes is proposed by using a more general mathematical model and a numerical 
technique based on the finite element method. The results presented here demonstrate 
that the model's behavior agrees with the behavior of biochemical reactions as it's con- 
sistent with the physical interpretation of the phenomena. Moreover, after comparison 
we can observe a complete consistency with literature findings [9,13-15,18]. A variety of 
numerical experiments were presented to confirm the accuracy, efficiency, and stability of 
the proposed method. In particular, the scheme was shown to be unconditionally stable 
and second-order accurate in space. 
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